{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Least Squares\n",
        "In this homework we find a method to find the least-squares solution using gradient descent with a constant step size. In this code, we will compare the closed form solution to the iteratively solved solution."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "#Run these to import all of the necessary files\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "import scipy.io as spio\n",
        "import scipy.signal as sig\n",
        "import scipy.io.wavfile as wf\n",
        "import matplotlib as plt\n",
        "from pylab import *\n",
        "import os\n",
        "import re\n",
        "import time\n",
        "from matplotlib import rc\n",
        "# Plots graphs in the browser window\n",
        "%matplotlib inline"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "#generate A matrix and y vector to use least squares regression\n",
        "#this code takes a while to run\n",
        "np.random.seed(121234)\n",
        "A = np.random.rand(30,2)\n",
        "y = np.random.rand(A.shape[0],1)\n"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## a)\n",
        "Create the function LS_closed_sol which takes in a matrix A and vector y and returns the square error using the closed form solution of least squares regression. \n",
        "As a reminder, the squared error is: $$||A\\vec{x}-\\vec{y}||^2$$\n",
        "The closed form solution for $\\vec{x}$ is: $$\\vec{x} = (A^TA)^{-1}A^T\\vec{y}$$"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def LS_closed_sol(A,y):\n",
        "    \"\"\"\n",
        "    The following function computes the least squares solution x_star\n",
        "    Inputs: \n",
        "    A : data matrix\n",
        "    y : desired output\n",
        "    \n",
        "    Ouputs: \n",
        "    e2 : least squared error possible for the given system\n",
        "    x_star : least squares solution for x\n",
        "    y_star : best approximation of y\n",
        "    \"\"\"\n",
        "    AT = A.T\n",
        "    ATA = A.T.dot(A)\n",
        "    x_star = np.linalg.inv(ATA).dot(AT).dot(y)\n",
        "    y_star = A.dot(x_star)\n",
        "    e2 = pow(norm(y_star-y),2)\n",
        "    return e2, x_star, y_star"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "LS_closed_sol(A,y)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## B)\n",
        "Given the step size $\\alpha$, we can create our function to iteratively solve for $\\vec{x}$.\n",
        "\n",
        "Create a function LS_iter(A,y,alpha,x2,iters) which takes as input a matrix A, vector y, constant alpha, vector x2, and constant iters. The function should calculate \n",
        "    $$\\vec{x}(t+1) = \\vec{x}(t)-\\alpha A^T(A\\vec{x}-\\vec{y})$$ \n",
        "up to time step $\\vec{x}(\\text{iters} - 1)$ as well as the squared error at each time step. The function should return an array which contains the squared error $||A\\vec{x}(t)-\\vec{y}||^2$, $x(t)$ and $y(t)$ for each time step."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def LS_iter(A,y,alpha,x2,e2,iters):\n",
        "    \"\"\"\n",
        "    The following function iteratively calculates x(t) for the \n",
        "    following discrete time system:\n",
        "        x2(t + 1) = x2(t) - alpha*A^T(Ax(t) - y)\n",
        "    input: \n",
        "        A : Data matrix\n",
        "        y : Desired output\n",
        "        alpha : step size in gradient descent\n",
        "        x2 : guess solution to equation A*x2 = y\n",
        "        e2 : the final squared error that we are going to get\n",
        "        iters : number of iterations for which x(t) is computed\n",
        "    output:\n",
        "        e : difference between squared error |A*x2 - y|^2 for each time step and the final e2 provided\n",
        "        x : x2(t) for all time steps\n",
        "        y : y(t) = A*x2(t) for all time steps\n",
        "    \"\"\"\n",
        "    at = A.T\n",
        "    ata = at.dot(A)\n",
        "    alata = alpha*ata\n",
        "    alaty = alpha*at.dot(y)\n",
        "    e = []\n",
        "    x = []\n",
        "    p = x2\n",
        "    estimates = []\n",
        "    for i in range(iters):\n",
        "        p = p-alata.dot(p)+alaty\n",
        "        e.append(pow(norm(A.dot(p)-y),2)-e2)\n",
        "        x.append(p)\n",
        "        estimates.append(A.dot(p))\n",
        "    return e, np.asarray(x), np.asarray(estimates)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## C) \n",
        "\n",
        "Now, let's plot the convergence of $x(t)$ and $y(t)$ as time progresses. Let's also plot their norms and the error as a function of time. \n",
        "\n",
        "Before you run the following code block, take a minute to guess how these plots would look. Do they make sense?"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "# First let's find the minimum and maximum eigenvalues\n",
        "l, v = np.linalg.eig(A.T.dot(A))\n",
        "l_min = np.min(l)\n",
        "l_max = np.max(l)\n",
        "\n",
        "# Next, let's calculate the least squares x_star, y_star and error e2\n",
        "e2, x_star, y_star = LS_closed_sol(A, y)\n",
        "\n",
        "\n",
        "\n",
        "# Finally, let's see how our solution x_iter and error evolves with \n",
        "# time/iterations\n",
        "alpha_arr = [1/l_min, 2/(l_min + l_max), 1/l_max]\n",
        "x2 = np.zeros((A.shape[1], 1))\n",
        "iters = 30\n",
        "\n",
        "for alpha in alpha_arr:\n",
        "    e_iter, x_iter, y_iter = LS_iter(A, y, alpha, x2, e2, iters)\n",
        "    \n",
        "    fig=plt.figure(figsize=(12, 4), dpi= 80, facecolor='w', edgecolor='k')\n",
        "    fig.suptitle('Gradient descent with step size = ' + str(alpha))\n",
        "    plt.subplot('131')\n",
        "    plt.title('Movement of $x$ towards $x^*$')\n",
        "    plt.plot(x_iter[:, 0], x_iter[:, 1], '.', color='tab:blue')\n",
        "    plt.plot(x_star[0], x_star[1], 'o', color='tab:red')\n",
        "    plt.xlabel('$x_1$'); plt.ylabel('$x_2$')\n",
        "    plt.legend(['$x$', '$x^*$'])\n",
        "\n",
        "    plt.subplot('132')\n",
        "    plt.title('Movement of $x$ with time')\n",
        "    plt.plot(np.arange(iters), x_iter[:, 0], '-', color='tab:orange')\n",
        "    plt.plot(np.arange(iters), x_iter[:, 1], '-', color='tab:purple')\n",
        "    plt.xlabel('time or iterations'); plt.ylabel('components of $x$')\n",
        "    plt.legend(['$x_1$', '$x_2$'])\n",
        "\n",
        "    plt.subplot('133')\n",
        "    plt.title('Gap in error $||\\epsilon||^2$')\n",
        "    plt.plot(np.arange(iters), e_iter)\n",
        "    plt.xlabel('time or iterations'); plt.ylabel('Error')\n",
        "    plt.yscale('log')\n",
        "    plt.legend(['$\\epsilon^2$ - closed-form-error'])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "scrolled": true
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Let's try this approach where the step size is even smaller! \n",
        "# and we can use more steps, to be near to a `continuous` approach to the solution\n",
        "alpha_arr = [.01]\n",
        "x2 = np.zeros((A.shape[1], 1))\n",
        "iters = 300\n",
        "\n",
        "for alpha in alpha_arr:\n",
        "    e_iter, x_iter, y_iter = LS_iter(A, y, alpha, x2, e2, iters)\n",
        "    \n",
        "    fig=plt.figure(figsize=(12, 4), dpi= 80, facecolor='w', edgecolor='k')\n",
        "    fig.suptitle('Gradient descent with step size = ' + str(alpha))\n",
        "    plt.subplot('131')\n",
        "    plt.title('Movement of $x$ towards $x^*$')\n",
        "    plt.plot(x_iter[:, 0], x_iter[:, 1], '.', color='tab:blue')\n",
        "    plt.plot(x_star[0], x_star[1], 'o', color='tab:red')\n",
        "    plt.xlabel('$x_1$'); plt.ylabel('$x_2$')\n",
        "    plt.legend(['$x$', '$x^*$'])\n",
        "\n",
        "    plt.subplot('132')\n",
        "    plt.title('Movement of $x$ with time')\n",
        "    plt.plot(np.arange(iters), x_iter[:, 0], '-', color='tab:orange')\n",
        "    plt.plot(np.arange(iters), x_iter[:, 1], '-', color='tab:purple')\n",
        "    plt.xlabel('time or iterations'); plt.ylabel('components of $x$')\n",
        "    plt.legend(['$x_1$', '$x_2$'])\n",
        "\n",
        "    plt.subplot('133')\n",
        "    plt.title('Gap in error $||\\epsilon||^2$')\n",
        "    plt.plot(np.arange(iters), e_iter)\n",
        "    plt.xlabel('time or iterations'); plt.ylabel('Error')\n",
        "    plt.yscale('log')\n",
        "    plt.legend(['$\\epsilon^2$ - closed-form-error'])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "scrolled": true
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.4"
    },
    "nteract": {
      "version": "0.23.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}